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Abstract The time dependence of deviations from the 
Gaussian state in a freely cooling homogeneous system of 
smooth inelastically colliding spheres is investigated by 
kinetic theory. We determine the full time dependence 
of the coefficients of an expansion around the Gaussian 
state in Generalized Laguerre polynomials. Approximat- 
ing this system of equations to sixth order, we find that the 
asymptotic state, where the mean energy T follows Haff 's 
law with time independent cooling rate, is reached within 
a few collisions per particle. Two-dimensional molecular 
dynamics simulations confirm our results and show expo- 
nential behavior in the high-energy tails. 



Keywords: Inelastic particles, velocity distribution, non- 
Gaussian behavior, high energy tails. 



Introduction 

Freely cooling systems of smooth inelastically colliding 
spheres or discs (in the following denoted by inelastic 
hard spheres systems (IHS)) have been investigated by 
means of kinetic theory and coniputer simulations by sev- 
eral groups (see e.g. |0J,i|,|,|,§). Most of the studies 
focus on latest times where interesting phenomena like 
formation of vortex patterns and clustering 
can be observed. For short times or not too high inelas- 
ticities, however, the system remains homogeneous with a 
decreasing temperature, or equivalently, a decreasing aver- 
age velocity. It is the so called Homogeneous Cooling State 
(HCS) the regime that will be studied in this paper. The 
HCS admits a scaling solution, i.e. if one scales all veloci- 
ties with the average velocity and assumes that the shape 
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of the scaled velocity distribution function remains con- 
stant in time, the entire time dependence is given by the 
time dependence of the average velocity only. This scaling 
solution is the starting point for a hydrodynamic analysis. 
Although many of the existing theories use a Gaussian ve- 
locity distribution function (which may be valid for small 
inelasticity) in general the shape is not Gaussian. First 
evidence, at very late evolution stages, was obtained by 
Goldhirsch et al. by measuring the fourth moment of 
the velocity distribution function. Later Goldshtein and 
Shapiro ^ proposed a solution based on an expansion in 
Sonine polynomials. Van Noije and Ernst correctly 
calculated the first term to linear order and Brilliantov 
and Poschcl pl included nonlinear corrections. Numerical 
solutions of the Boltzmann equation were calculated us- 
ing Direct Simulation Monte Carlo method (DSMC) |2|. 
Finally, extensions to viscoelastic particles have been re- 
cently presented in jl^ . 

Further confirmation of non Gaussian behavior is given 
by Esipov and Poschel by studying the high-velocity 
tails of the velocity distribution function. They find that 
the tails are of an exponential type instead of a Gaus- 
sian, results confirmed by DSMC by Brey et al. and 
also in experiments performed by Losert and coworkers 
[Em. In fact, an asymptotic exponential tail of the form 



exp(- 



with (3 < 2 seems to be a more fundamental 



behavior, as it is also present in driven or vibrated granu- 
lar materials Recently, a detailed study by DSMC 
in driven systems with three different types of forcing has 
been presented |T^. Only for the so-called Non-Gaussian 
thermostat (where there is a balance between energy in- 
put and dissipation) [3 = 2, while for the other forcings, 
/? = 3/2 and 1. Surprisingly, no molecular dynamics re- 
sults have been presented so far for calculating moments 
and high energy tails for the freely evolving case. 

Our starting point is similar to i.e. expand the 
scaled velocity distribution function in a series of General- 
ized or Associated Laguerre polynomials around the Gaus- 
sian distribution with coefficients denoted with a;. How- 
ever, we assume that the coefficients are time dependent 
[ pdl . With these ideas we try to achieve two goals. Firstly, 
investigate the influence of higher coefficients as, . . . , ae, 
in the expansion of the distribution function. Secondly 
study their time evolution. We find that for not too high 
inelasticities the above expansion seems to be convergent. 
Furthermore, the cooling proceeds in two stages: (1) A 
fast decay (in the order of few collisions per particle) of 
all coefficients a; to their asymptotically constant values. 
(2) An algebraically slow decay of the kinetic energy T 
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determined by Haff's law -^T — CT^'^^'^ and time inde- 
pendent coefficient C depending on the asymptotic values 
of a;. Two dimensional event driven molecular dynamics 
simulations are performed in order to test the theory with 
good agreement for moderate inelasticity. For higher in- 
elasticity the perturbations expansions seems to fail and 
in the simulations we were able to observe a transition to 
an exponential high-energy tail. 

The paper is organized as follows: In Sec. 2 we pro- 
pose the expansion of the velocity distribution in Laguerre 
polynomials with time dependent coefficients. In Sec. 3 we 
determine formally the full time dependence of the HCS 
expressed by the time evolution of the coefficients of the 
expansion. We obtain an infinitely large system of ordi- 
nary differential equations, which can only be investigated 
approximately. This is done in Sec. 4, where a truncation 
scheme is proposed and analyzed to different orders. In 
Sec. I we compare the analytical theory to results from 
event-driven simulations and the validity of the pertur- 
bation expansion is discussed. Results for the exponential 
high-energy tail are presented here. We summarize the re- 
sults in Sec. |[ 



The system under consideration 

We consider a system of N smooth, inelastically collid- 
ing spheres with diameter a confined to a d-dimensional 
volume V , so that the homogeneous density is given by 
n := Y- The positions of each sphere are denoted by 
and each particle has a velocity Vi. The particles inter- 
act via a hard-core potential and in each collision (i.e. if 
rij :— \rij\ := jr^ — rj\ = cr) the velocities are instan- 
taneously changed by the following collision rules deter- 
mined by a constant coefficient e„ € [0, 1] of restitution 



2 

1 + e„ 



{vij ■ rij)rij , 



(1) 



where Vi 



Vj and rij 



rij /rij. Velocities after 



collision are primed quantities given by velocities before 
collision (unprimed quantities). 

The IHS is described statistically by the single particle 
distribution function p(r,v,t)drdv, the (average) number 
of particles at positions between r and r + dr and with 
velocities between v and v + dv at time t. As proposed 
in 1^ , for a homogeneous inelastic system the distribution 
function can be expressed by a scaling function as: 



1 



pic,t), 



(2) 



where c — v /vq (t) and vq (t) is the thermal velocity defined 
as the square root of the second moment of the distribu- 
tion function: 



dvp{v,t)v'^ = n-vl{t) . 



(3) 



The temperature is then defined as T{t) :— y'^oW- For 
elastic systems the distribution function is Gaussian and 



it is expected that it will remain close to a Gaussian for 
small inelasticity. Therefore, we expand the scaled distri- 
bution function in a series of Generalized or Associated 
Laguerre polynomials around the Gaussian distribu- 
tion function. The expansion is carried out in the scaled 
velocity variable c and with time dependent coefficients 
ai (t) . The general ansatz for the single particle distri- 
butions function for the homogeneous cooling then reads 

oo 

p(c, t) := exp {-c')J2 (^)^r (c') , (4) 

1=0 

where a = d/2—1 in d dimensions. In the context of kinetic 
theory Laguerre polynomials are called Sonine polynomi- 
als II . 

The normalization condition for p, J dvp = n, leads to 
ao = 1. We express by the first and second Laguerre 
polynomial 



(5) 



and using the orthogonality relations for the Laguerre 
polynomials we find 



J dv p{v, t)v^ = '^"o " I "1 



(6) 



which iniplics together with Eq. (^) that ai — for all 
times We denoted the binomial coefficients by 

it)- . 

Finally, as Laguerre polynomials are orthogonal, the 
coefficients ai are given by 



nai{t) 



1 



cr) 



dvp{v,t)L1 



Mi). 



(7) 



The Homogeneous Cooling State 
The Boltzmann Equation 

We assume that dynamics of the one particle distribution 
function p is given by the Enskog Boltzmann equation, 
which can be written in d dimensions without external 
forces as 

dtp{r, vi,t) + {vi ■ Vr)p{r, f i, t) = J[p, p] , (8) 
with collision integral 

J[/9, p] = cr'^^^x j dv2 j d&0{vi2 ■ ar){vi2 ■ or) 



-^-l]{pir,Vi,t)p{r,V2,t)). (9) 



& is the unit vector pointing from particle 2 to particle 
1, X the pair correlation function at contact, and de- 
scribes 'restituting collisions' by changing velocities in p, 
i.e. b~^p{r,v" ,t) = p{r,v,t) in a way that v" are the 
velocities before collision leading to v after collision. The 
operator b describes 'direct collisions' given in Eqs. (0). 
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The inverse operator of b, i.e. b~^, is simply given by sub- 
stituting e„ by l/e„ in Eqs. (0). 

By multiplying the Boltzmann equation Eq. (H) with 
some function ijj{vi) and integrating over Vi one gets 

dt J dviilj{vi)p{r,Vi,t)+V- J dvivi'ip{vi)p{r,vi,t) = 

dv,4,{v,)J[p,p] , (10) 



which can be rewritten in the form of a balance equation 

5tV^+v-jv = 9rV, (11) 

describing that the time change of an averaged quantity 
ip is due to flux j'^ or due to change through collisions. 
The right hand side of Eq. (|l0| ) can be written as [|l8| 

J dvMvi)J[f, f]^a''-\ J dv^dvi J da 

e{vi2 ■ ir){vi2 ■ a)p{vi)p{v2)A^ , (12) 

and Aif) is the change of ^ in a direct collision for both 
particles Aij; — \{'ij){v\) + iIj{v2) — — ^'(■"2))- 



Dynamics of Moments 

Using Eqs. (||) and ( p^ the time dependence of T{t) 



in the homogeneous case is given by 



-T 



d m 



dt dt 2 
where 7 is defined as 



-27^0^ , 



(13) 



/27r 1 



7 



dSd 



-ij j dcidc2d&0{ci2 ■ o-)(ci2 • it) X 

p(ci)p(c2)(6-l)i(c? + c2) , 



(14) 



and p as in Eq. (^. Sd is the surface of a unit sphere in 
d dimensions and ujq the Enskog collision frequency for a 
classical gas of hard spheres with temperature T, given 
by: 



r{d/2) 



and loq 



Sd 



Xn{2ar-'vo . (15) 



If the velocity distribution function p in Eq. (Q) is a 
Maxwellian, 7 takes the value of 70 := (1 — ef^)/{2d). This 
is the Gaussian value of the energy decay rate obtained by 
Haff |l^ . The fact that this distribution function is not a 
Gaussian modifies the cooling rate, as it will be calculated 
later on. 

In order to obtain the time evolution of a; , we take the 
time derivative of Eq. (^ , and it has to be considered the 
time dependence of p{v, t) as well as the time dependence 
of ^r((;;;7(t) )^) ■^i^ i'o(i)- The time dependence of p{v,t) 
is given by Boltzmann equation, and the time dependence 



of vo{t) is given by equation (|13[). After a straight for- 
ward calculation using differential formulas for the La- 
guerre polynomials, we get 



d 

dt"' 



uoli + l2-fujo{ai - a;_i) , 



(16) 



and 
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'2t: 1 1 



Sd 



1 /" 

— / dCi_dC2d(Te{cx2 ' 't)(ci2 • (t) X 



p(ci)p(c2)(6-l)-(Lr(cf) + Lr(c^)). (17) 



All collision integrals 7 and 7; depend on ai for all I via 
p. We mention here that our approach is equivalent to 
the dynamics proposed in jl^, but has the advantage to 
give immediately the explicit time dependence of all coef- 
ficients, at least formally. 



Collision frequency 

The set of equations (|l^) and ( [l6| ) for T and a; with 7; and 
7 given by (^) and ( |14[ ) are the main results of this pa- 
per. However, before analyzing them in detail in the next 
section and comparing them with computer simulations 
in Sec. 5, it is instructive to study the collision frequency 
ijj. In order to do so, we introduce the average number 
of collisions, r, that a particle has suffered in a time t. 
Then, the collision frequency is defined as w = -^T{t). In 
elastic fluids a; is a constant number depending only on 
the density and temperature, so that r and t are propor- 
tional quantities. In granular fluids, however, oj depends 
on time, as the temperature (and more precisely also the 
shape of the distribution function) of the system changes 
with time. Therefore, it is more natural from a physical 
point of view to express the time evolution equations in 
terms of the variable r rather than t. Moreover, the hy- 
drodynamic matrix become time independent when the 
l^drodynamic equations are expressed in the variable r 

To determine uj = jt'''{t) we use Eq. ( [T^ ) and the fact 
that in each collision the number of collisions that each 
particle has suffered increases by one and we obtain 



d_ 
It 



T = ujoIt , and 



(18) 



It = 



/27r 1 

's7^ 



dcidc2da-0{ci2 ■ o-)(ci2 • (t)/5(ci)/5(c2) . 

(19) 



7r depends on all a; and for the case that all a; = for 
Z > 1 we would get 7r = 1 and thus the Enskog value loq. 
We deflne a time f by 



df = LOodt . 



(20) 



Note that f is only an approximation of r defined in 
Eq. (|l8|), so it does not really measure time in collisions, 
but we will show later that the deviations of f from r re- 
main small for not too high inelasticities. In other words 
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we hope that the colhsion frequency is approximately de- 
termined by the Enskog value and corrections due to devi- 
ations from the Gaussian affecting the collision frequency 
are small. 



Cooling rate 

How can the dynamics be described in a state where all 
coefficients have already reached their asymptotic values? 
Note that the quantities 7 and 7,- are entirely given by 
the values of a; . Assuming that all a/ have reached their 
asymptotic values for some time t > t* or equivalently 
T > T* the quantities 7 and also remain constant and 
we denote their asymptotic values by 7* and 7* . Then we 
consider Eq. (|T|) and (|l|): 



dt' 



-T 



-2-f*LOo{T)T , 



dt 



(21) 
(22) 



which is solved analytically 



T 



T{t* 



T(t*) exp(-2777;(T 



r*)) , (23) 



so that 



r{t) 



7* 



ln[l + 7*c.o(r(r))(t-i*)] . (24) 



Eq. (|24| ) provides a relation between coUisional time and 
real time. Eq. (p3|), i.e. the algebraic decay of the temper- 
ature in time or the exponential behavior in r is called 
Haff 's law. Furthermore, this equation makes explicit the 
fact that the shape of the velocity distribution function 
modifies the energy decay rate. 



We will concentrate on two aspects: (i) To investigate 
the dynamics we integrate the full set of differential equa- 
tions (up to a certain order A). The asymptotic values of 
the coefficients can then be obtained by taking the long- 
time limit if they become stationary in time, (ii) To discuss 
the stationary state we set the left hand side of Eq. ( p!^ ) 
equal to zero. This set of coupled and, as the case may be, 
non-linear equations can be solved with the numerical tool 
provided by the computer algebra program. Note that not 
all of these stationary values are necessarily dynamically 
stable solutions of the corresponding differential equation. 



Results to order 2 

In a first step we only take into account 02 to linear order. 
Then the functional form of the equation for 02 Eq. ([l^ ) 
is given by 



^a2=a;o(72 



- 4702) -7^02 = 72 

dr 



4702 



(25) 



Again, the use of f simplifies the form of the equations, 
as it eliminates the time dependent factor uq. We recall 
here that both 72 and 7 depend on all a;, but for the 
approximation treated here, they only depend on 02 in a 
linear manner. Therefore, we can express Eq. (^) as 



d 



A + Ba2 + 0{\'), 



(26) 



where A and B are constants given by the collision in- 
tegrals in 7 and 72 with the explicit expressions in two 
dimensions: A = |(e^ - l)(2e^ - 1) and B = ^^(SOe^ - 
5e2 -32e„-57). 

(i) Dynamics- The evolution equation ( p6| ) for 02 where 
time is expressed in collisions per particle is linear in 02, 
so it can be easily integrated to give, when 02(0) = 0, 



Analytical Results 



a2(T) = -^(l-cxp(i3f)) , 



(27) 



Truncation scheme 

Up to now we have determined the full time dependence 
of the HCS in terms of the time dependence of all its mo- 
ments in Eqs. ( p^ ) and (|l^). This infinitely large system 
of differential equations can only be solved by truncation. 
The approximate solution found by truncation only makes 
sense if all neglected terms are small as compared with the 
remaining ones. On the other hand, for small inelasticities 
the velocity distribution function is close to a Gaussian, 
so that 02 <C 1 = flo. We generalize this unequality and 
assume that a;+i <C a/ for all i.e. contributions from 
higher order coefficients get smaller the higher the index. 
Therefore, we propose a truncation scheme in which we as- 
sume that ai is of order A', where A is a small parameter. 
If we now make "an approximation of order ©(A')" we ne- 
glect all terms in all considered equations higher than A' . 
This truncation scheme produces a finite set of differential 
equations that can be solved. 



so that the asymptotic value of 02 is reached exponentially 
fast on a time scale of the order of fg := —B~^. The decay 
time To ranges between 1.7 and 2.25. As an important 
consequence the asymptotic solution is quickly reached on 
a kinetic time scale of few collisions per particle. 
(ii) Stationary state- For times larger than tq 02 reaches 
the stationary value of —A/B which coincides with the 
values calculated in iflCll. 



Results to order 3 

To keep the discussion simple and to compare results from 
order to order, we first take into account only 02 and 03 
i.e. up to ©(A^). We then still have to deal with equations 
which are linear in the coefficients (a| is already of order 
0{X'^)). In the next section we will discuss the non-linear 
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case up to order O(A^). The equations read 
d 



dt 
d 



dt 
d_ 
dt 



wo72 + 47CJo(a2 



0.2 



as = wo73 + 671^0(03 



0), 
02) 



dt' 
d_ 
'dt 



(28) 



f = cjq neglecting corrections of a2 and 03. 



We use computer algebraic programs to calculate the col- 
lisions integrals 7, 7; and 7r up to order 0{\^). The ana- 
lytical solutions are rather lengthy and we will only show 
here results for a system with e„ = 0.9 in Figs. 1-3. 
(i) Dynamics- We have solved the simultaneous Eqs. ( p8| ) 
numerically^ for the case e„ — 0.9 and in the following we 
always plot time in units of 1/ujq{T{Q)) and temperature 
in units of T(0). We have chosen 02(6) = a3(0) = as 
initial condition. In a first step we proof that the approx- 
imation to use f instead of t can be justified (at least to 
this order). In Fig. || a) we show the relative deviations of 
the true number of collisions to the approximation given 
by the Enskog Boltzmann value, i.e (r — f)/T as function 
of time in a semilogarithmic plot. 
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Fig. 1. a) Relative deviations of collisions per particle from the 
approximation given by the Enskog value as a function of time, 
b) Collisions per particle as a function of time. Dissipation, 
e„ = 0.9. 



We see that the relative deviations remain smaller than 
0.2 %. This allows us, at least in the homogeneous cooling 
state, to use f instead of r in Eqs. ( ^ and (^J). 

In the asymptotic state we get 7* — 0.04723 for e„ — 
0.9 and values from the numerical integration of Eqs. ( p8| ) 
coincide with Eq. ( |2^ ) and (24) within the graphical accu- 
racy, so we plotted here only the numerical solution. Fig. 
|l| b) shows rit) which has the same form as predicted in 
Eq. H). 

In Fig. a), we show T as a function of time in a double 
logarithmic plot. We see the well-known asymptotic time 
dependence T cx t^"^. In Fig. || b), we show T as a function 
of T in a semi logarithmic plot resulting in a straight line 
with slope —27* as predicted by Eq. (|2^). 

In Fig. ^ we show the time dependence of 02 and 03 as 
a function of time a) and as a function of r b). 

^ We have used the built-in numerical procedure dsolve of 
MAPLE to integrate the differential equation. 
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Fig. 2. Temperature as a function of time a) and collisions per 
particle b). Dissipation, e„ = 0.9. 
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Fig. 3. Coefficients 02 and a-i as a function of time a) and 
collisions per particle b). Dissipation, en = 0.9. 



We see that 02 and 03 reach their asymptotic value 
on a very short time scale which is of the order of few 
t's. Therefore, few collisions per particle are necessary to 
reach the asymptotic state for 02 and 03. 
(ii) Stationary state- As mentioned above we calculate 
the stationary values by setting the l.h.s of Eq. (|2^) equal 
to zero. In Fig. |^ we show the results for the stationary 
values of 02 and 03 to 0{\'^) as well as 02 to 0{}?). 
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3.2 to order 3 
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^ ^ 85 to order 2 
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Fig. 4. Coefficients 02 to 0{\^) 
function of e-n 



and 02 and 03 to 0{\^) as a 



As long as e„ > 0.6, 02 to 0{\^) does not differ sig- 
nificantly from 02 to 0{\^) and 03 remains small. We see 
stronger differences for smaller e„ and 03 becomes as im- 
portant as 02 which indicates stronger deviations from the 
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Gaussian state. We also cannot assume anymore that cor- 
rections of higher orders remain small since we do not 
have any indication that the series is converging in the 
sense that the |a;| are small and decreasing. 

Since to 0{X'^) we have to deal with a set of linear 
equations we only find one unique solution. Considering 
higher orders one will find many solutions whose validity 
must be investigated. We will discuss this problem in the 
next section. 



Results to order 6 

In this section we go to 0{X^), which is the highest order 
we were able to calculate with the computer algebra pro- 
gram. 

(i) Dynamics- In Fig. |^ we show for e„ = 0.8 the dy- 
namics for the 5 non- vanishing coefficients a2, ■ ■ ■ , og as 
a function of time. We have chosen the initial condition 
02(0) = ... = ae{0) = 0. We see again a very fast decay 
to their asymptotic values. We observe that \ai\ > |a;+i| 
and in this sense the perturbation expansion seems to con- 
verge. 
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Fig. 5. The coefficients 02,..., ag calculated to order 0{\^') 
as a function of time for e„ — 0.8. 



(ii) Stationary state- We calculate the stationary val- 
ues by setting the l.h.s of Eq. (|l6|) equal to zero. In Fig. || 
we show the results of the stationary values as a function 
of e„ > 0.3. 

For e„ > 0.7 the coefficients remain small and the 
expansion seems to converge in the sense that |a;| > |a;+i| 
for all I. For e„ < 0.7 the absolute values of the coefficients 
start to grow and seem to diverge with e„ approaching 0.3. 

To discuss the validity of these results we compare in 
Fig. 1^ the asymptotic values of 02 of order 0{\^) up to 
order 0{X^). As long as e„ > 0.6, we do not find significant 
differences between the two orders, only the first order 
differ slightly from the other ones. This is a further hint 
that for these values of e„ the perturbation method works. 
Moreover, the ratios of a/^i/a/ are small and of the same 



0.08 




_0.04 I ' ' ' ' ' ' ' 1 

0.6 0.7 0.8 0.9 1 

Fig. 6. Stationary values 02, . . . , as calculated to order 0{X^) 
as a function of e„ 

order: for instance, for e„ — 0.85, these ratios are: 03/02 — 
0.39, 04/03 ^ 0.29, 05/04 = 0.24 and 05/05 = 0.16. 

For e„ < 0.6 the results differ drastically from order 
to order and the proposed truncation scheme for Eq. ( [l^ ) 
fails. We conjecture that around ~ 0.6 an essential 
change in the distribution function occurs. Then the dis- 
tribution function is no more described by small deviations 
around a Gaussian and might be better expressed by an 
expansion around an exponential as suggested in and 
confirmed by DSMC simulations of iQ. We will go back 
to this point at the end of Sec. 5. 
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e„ 

Fig. 7. Stationary value of 02 to order 0{\"*) up to O(A^) as 
function of Cn- 



In view of Fig. |7| another possible explanation for the 
failure of the convergence is that the expansion in A is 
of asymptotic type, in such a way that including higher 
orders the expansion would break down in the whole range 
< e„ < 1. Unfortunately, we cannot decide which is the 
correct option, as we can only calculate up to order 0{X^). 
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Further unstable solutions 

Since we have to deal with non-linear equations, the so- 
lution is not unique and e.g. for e„ — 0.8 two further 
stationary solutions can be found, similar as in We 
list the values of the other coefficients for these stationary 
solutions: 





Solution 1 


Solution 2 




8.95 


-24.62 


as 


-f4.39 


-4.50 




59.11 


39.53 


as 


-109.17 


178.6 




127.8 


-197.7 



Both solutions are dynamically unstable which we have 
shown by numerical integration of the corresponding dif- 
ferential equations ( [l^ ) and (|l^). In addition we observe 
that the higher coefficients are not at all negligible so 
that our assumptions, which should allow us to truncate 
the system of differential equations, are severely violated. 
Hence for these cases we can assume that we have not 
even found an approximate solution of the homogeneous 
Boltzmann equation. 



Computer simulation results 

In the literature only Direct Simulation Monte Carlo meth- 
ods (DSMC) have been used for measuring the values of 
02 and 03 |0 and 02 agrees very well with the value cal- 
culated in However, DSMC lacks some features of the 
real IHS fluid, as correlations among the particles. 

We present in this paper for the first time results for 
02 and for high-energy tails obtained from Molecular Dy- 
namics (MD) simulations of the IHS system in 2 dimen- 
sions. Our code closely follows the event driven molecular 
dynamics code presented in ||2^, adapted to the collision 
rules described in Eqs. (Q) and accelerated by techniques 
described in [ pT| . Typical simulations are performed with 
N ~ 50000 particles in a square box of size L, being its 
area fraction (j) = — ^^^f . The initial configura- 



tion is that of an elastic fluid at equilibrium (Maxwellian 
distribution for velocities and equilibrium correlations for 
the positions due to excluded volume effects) , prepared by 
running the system with e„ = 1 (elastic interactions) for 
not less than 50 collisions per particle. Therefore, in the 
initial state o; = for Z > 1. 

At the beginning of the inelastic evolution the system 
remains for some time in the HCS, where the assumptions 
made in Sec. 2 are fully applicable and where computer 
simulations will serve to test those predictions. Later, vor- 
tices and clusters start to develop through the system and 
homogeneity is lost ]l|,D. The higher the density cj) and 
the inelasticity the sooner these structures appear and 
important deviations from the theoretical values of ai are 
expected. We will come back to this point later. Further- 
more, the analytical results are independent of the density^ 

^ To eliminate the dependency of real time on the density, 
time can be scaled by the collision frequency uj{T{0)) at time 
t = 



when expressed in r, but only depend on the inelasticity 
e„. Hence we have performed our simulations at low den- 
sity of (j) = 0.05, although simulations at higher and lower 
densities have also been carried out. 



Results for moments 

The typical time evolution of the 4th cumulant is shown 
in Fig. ^, where we have plotted the value of 02 versus the 
number of collisions per particle r for a low density case 
(f) = 0.05 and low inelasticity e„ = 0.92. The dotted line 
is the result of Eq. ( p7| ) while the solid line is the result 
of the numerical simulation averaged over two realizations 
to slightly improve the accuracy. In this plot we observe 
the typical features of the IHS evolution described in for- 
mer sections. Initially 02 is equal to zero, as the system 
starts from a Maxwellian distribution, with o; = for 
I > 1. Then, within a very short time of a few collisions 
per particle, deviations from a Maxwellian build up in the 
system and the asymptotic values of o; are reached. This 
is a very fast process on a hydrodynamic time scale, as it 
involves only a few collisions per particle and, therefore, 
a few mean free times. Then the moments stay constant 
(within the accuracy of the computer simulations) as long 
as the system remains in the HCS. 



-0.01 



-0.02 - 



-0.03 



-0.04 




Fig. 8. Time evolution of a2 for an IHS system with (p — 0.05 
and e„ — 0.92. The dotted line is the analytical result, while 
the solid line is the numerical simulation data averaged over 
two realizations. 



A best fit of the simulation data to an expression like 
Eq. (13) gives that = 1-6 ± 0.5 and 02(00) = -0.026 ± 
0.004, while the theoretical values developed in Sec. 4 are 
To = 1.85 and 02(00) = —0.0246. We observe excellent 
agreement with the theory. Unfortunately, the accuracy of 
our computer simulations is not high enough to distinguish 
between the lowest order and the order 6. The DSMC 
method also finds good agreement with the value of 02 

At later times the system is no longer in the HCS 
and the assumptions used in the theoretical sections break 
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Fig. 9. Time evolution of 02 for an IHS system with (j> = 0.03 
and fin = 0.4. The dotted line is the analytical result, while the 
solid line is the numerical simulation data averaged over two 
realizations. Deviations from HCS are seen after r > 12. 



Fig. 10. Coefficient 02 versus the coefficient of restitution 
e„. The solid line is the theoretical prediction of Eq. ( |26| ) and 
the circles are the values calculated from MD simulations with 
their corresponding error bars. 



down. This is best illustrated in the Fig. y, where a simu- 
lation at low density — 0.03 but at very high inelasticity 
Cn = 0.4 is presented. We observe the same features de- 
scribed in Fig. ^ with values from simulations tq — 1.5±0.5 
and 02(00) = 0.125±0.007, a best fit for r < 10, compared 
to the theoretical values tq = 1.83 and 02(00) = 0.130. 
Again, the agreement is excellent. However, after a short 
time T > 12 the values of the moments start to deviate 
from the theoretical predictions. For r > 12 the homo- 
geneity assumption breaks down. This can be checked, 
e.g. by plotting the curve of energy vs time (2[||]: devi- 
ations from Haff's law imply lack of homogeneity 
Visual inspection of the system (not presented here) show 
that, at this high inelasticity, the system immediately de- 
velops currents and dense clusters where particles move 
almost parallel inside them. The description in terms of 
p(c) is wrong, as it does not take into account the local 
macroscopic currents. The values of 02 can, at this late 
evolution stages, grow up to values of 02 = 2 Q. This 
regime is outside the scope of this article. 

Concerning the duration of HCS, a hydrodynamic anal- 
ysis 1^ shows that currents develop through the system on 
a time scale of r of the order of 7q"^. In Fig. || the devia- 
tion with respect to 02(00) appears at r ~ 35 — 1.37(^^, 
while in Fig. ^ it appears at r ~ 12 = 2.h^^^. 

The results of the MD simulations for 02 in the HCS 
compared with the asymptotic solution of Eq. ( p6[ ) are 
given in Fig. 10. The agreement is excellent even down to 
high inelasticities as e„ = 0.4, as shown in Fig. y. However 
our simulations do not have precision enough to test if 
higher order corrections are important. Concerning the 
values of fo, MD results are very close, but always smaller 
than the theoretical values quoted below Eq. (^ and they 
are affected by large errors. 

To conclude this section, we have also measured the 
6th moment, (c®) = {v^) / {v'^)^, and from this quantity we 
have obtained 03 given by 03 = — ^(c^) + (1 + + 



\48 



7d 
12 



1). The results also agree with the solutions 



of Eq. (p8|), but, as expected, the discrepancies are now 
larger because the absolute value of 03 is very small. 



Cooling rate 

Another way to test the analytical results of Sec. 4 is 
to measure the dissipation rate given by the decay of 
the temperature T or energy E versus r. Haff's calcu- 
lations predicts an exponential law exp(— 27ot) with 70 = 
(1 — e\)/2d while higher order corrections modify it as 
shown in Eq. (|2^). These corrections respect to 70 are 
very small of only few parts in a thousand. If the IHS sys- 
tem is too large or inelastic, these deviations cannot be 
measured, as the curve of energy vs r bends apart of the 
exponential law (|,^. This is due to the appearance of 
currents and vortices as quantitatively explained by [p2| . 
Moreover, temperature and energy are no longer propor- 
tional, and, in contrast to energy, temperature is difficult 
to measure in a numerical experiment. 

However, if the system is small enough no shear or clus- 
tering instability is excited (see e.g. Refs.||,0| for detailed 
explanations) and the system is forced to remain in the 
HCS for all times, where Eq. ( p^ ) is valid. This is called 
'kinetic regime' in Ref.[|j. The drawback of this method 
is that, for a given density, it sets a maximum number 
of particles, that decreases with increasing inelasticity. As 
described in Ref. a lower bound for kmin = x ^1 
to be satisfied, in the notation of 0. We keep kmin = 2k\. 
This condition, together with our chosen low density of 
(j) — 0.05, restricts the number of particles to = 80 at 
e„ = 0.70, and the results are no longer reliable. Even 
for e„ — 0.95 the maximum number of particles is only 
N ~ 300. Hence we restrict ourselves to e„ > 0.7. 

Following this method we have performed simulations 
and have measured the energy decay rate as a function 
of T. We have verified that it is indeed exponential for all 
times, and, in order to improve the statistics, results are 
averaged over 1000 realizations. The results are presented 
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Fig. 11. Deviations of tlie cooling rate as a function of tlie 
inelasticity as obtained in simulations of small systems that 
remain in the HCS. 



in Fig. |ri|, where we plot the difference between the sta- 
tionary energy decay rate 7*/7^ from Eq. ( ^3| ) and 70. 
Circles are data obtained by MD simulations with their 
errors bars, while dotted and dashed lines are the results 
from our theory to order 0{X^) and 0{X'^) respectively. 
The small number of particles does not allow to obtain 
any significant result below e„ < 0.7. For e„ > 0.75 we 
observe a reasonable agreement, although we find signif- 
icant deviations with respect to the theoretical results, 
MD data are always smaller that theoretical values. Un- 
fortunately, no direct measurements with DSMC of the 
deviations respect to 70 have been reported so far. They 
would allow us to compare our results and elucidate the 
nature of these deviations. 

It is important to note here that due to the small size 
of the deviations of the cooling rate with respect to 70 it is 
necessary to use in the theory the real number of collisions 
T instead of f . The approximation f ~ r is correct within 
a few parts in a thousand as shown in Sec. 4, which is of 
the same order of magnitude as the correction j* /j^ with 
respect to 70. This is not the case, however, for calcula- 
tions of the dynamics of a/, for example given in Eq. (p7|), 
where we are not interested in such small deviations but 
want to give a first estimate of the time scales. Hence the 
use of f and the approximation made in Eq. (27) is fully 
justified. 

Finally, there is still the open question if deviations 
from the theoretical 02 and the cooling rate are due to the 
existence of correlation in the HCS and the breakdown of 
the molecular chaos hypothesis in Eq. (^ reported in 

These effects cannot be tested in DSMC either, as this 
method is based on the factorization of the two particle 
distribution function. 



High-energy tails 

In Rcfs. it has been shown that strong deviations 

with respect to the Maxwellian are present in the tails of 
the distribution, where particles have large energies. More 
precisely, they have found if c = u/wo(i) ^ (1 — e^)~^ the 
velocity distribution function is no longer Maxwellian, but 
a simple exponential p ~ ^exp(— Ac) instead. This expo- 
nential distribution has been verified by numerical solu- 
tions of the Boltzmann equation using the DSMC method 
[0. The shape of the exponential is very different to 
the Maxwellian and therefore it is understandable that 
an expansion of the type given in Eq. (^) might be non- 
convergent as suggested by our theoretical analysis. An- 
other reasonable possibility is that the series is indeed 
convergent but we have not gone high enough in the trun- 
cation scheme or there is a better choice of truncation, 
because a;+i is as large as ai, as shown in Fig. 6. 

In order to investigate the velocity distribution func- 
tion and make the exponential range accessible, we have 
performed MD simulations with extreme inelasticities of 
e„ = 0.1, 0.2, 0.4 to compare to moderate inelasticities 0.6 
and 0.8. Moreover, as we need high accuracy in the tails, 
where populations are small, we have simulated systems 
with 250 000 particles at (p = 0.05. However, even with 
this large number of particles we are not able to obtain 
the accuracy that can be achieved by the DSMC method 
. We will only be able to give evidences of the exponen- 
tial tail. We measured the distribution function at times, 
where 02 has already reached its asymptotic value, but 
the homogeneity assumption is still valid. In the example 
of Fig. H this would be at times 5 < r < 10. 

To estimate the distribution function from data, we 
use the kernel estimator technique described e. g. in [ p5[ . 
In general, given a set of outcomes Xi, j = 1, . . . , M of a 
random experiment the distribution function p{x) can be 
estimated by 



1 



M 
i=l 



1 



1^8 



exp 



2,52 



(29) 



The idea behind this method is that each data point also 
gives some information about its surrounding, which can 
be justified for smooth distribution functions. It is not 
necessary to use a Gaussian as kernel in Eq. (|2^), any 
normalized and more or less sharply peaked function can 
be used. The value i5 is a free parameter which was chosen 
such that the measured distribution function for the elas- 
tic gas (i.e. the initial condition) is fitted best (to the eye) 
to the Maxwellian. We choose 5 = 0.05. Since the distri- 
bution function p((?^ -with c the modulo of the velocity- 
is not continuous at c = 0, i.e. p = for c < and 
p(0) « I/tt this technique gives bad results around c = 0, 
but better results than the histogram method for the in- 
teresting high- velocity limit, where only few data points 
are given. 

The simulation results for the velocity distribution func- 
tion p as a function of c as well as the Maxwellian are 
shown in Fig. 12 in a semi-logarithmic plot. For c > 3.5 



the statistical accuracy is poor, so results are only plotted 
up to c = 3.5. 
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Fig. 12. Mcixwellian and measured velocity distribution func- 
tion p as a function of c for various values of Cn. 



We observe that the deviations with respect to the 
Maxwellian are larger for lower values of e^. On the con- 
trary, as e„ increases, the measured distribution approaches 
the Maxwellian. Closer inspection shows that for e„ — 
0.1 the distribution gets possibly close to an exponential 
(straight line in the semilogarithmic plot of Fig. |2|) for 
2 < c < 3.5, while for e„ — 0.6 the range where logp 
seems to be linear shrinks to 3 < c < 3.5. We will show 
below that in this case (e„ = 0.6 and c < 3.5) the distri- 
bution function can be reasonably well described by the 
results of the perturbation expansion around the Gaussian 
given in Eq. . If we perform a linear fit to an exponential 
in these ranges, we obtain values of the coefficient A quite 
close to those reported by 0, tested in by DSMC 
method. For instance, for e„ = 0.1, A ~ —3.2, that in- 
creases to A ~ —3.8 at e„ = 0.4 and further to A ~ —4.7 
at e„ = 0.6. 

If we go beyond e„ > 0.6 the perturbation expansions 
of Sec. 4 seems to converge and it make sense to compare 
the analytical results with the simulation data. In Fig. ^ 
we show for e„ = 0.6 and e„ = 0.8 results of the simula- 
tions and of the analytical theory to orders and A® . 

In the whole range of c which is plotted, the analytical 
results to order A^ coincides fairly well with the measured 
data. However, results to order A^ agree at small veloc- 
ity, c < 3 for e„ — 0.8, but fails for higher velocities. It 
seems reasonable that higher orders are needed to describe 
higher energy tails, where deviations respect the Gaussian 
are larger. 

Therefore, for e„ > 0.6 there is no need to describe 
data by an exponential tail for high velocities, although it 
cannot be assured, that the theory does not fail for even 
higher velocities. Note that the distribution function at 
c « 3.5 is already smaller than 10"^ so that for 250 000 
particles only 2 or 3 particles might not be correctly de- 
scribed. 

6 

Conclusion 

In this article we investigated by means of analytical the- 
ory and simulations the dynamics of a freely cooling sys- 
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Fig. 13. Measured velocity distribution function p as a func- 
tion of c for Cn = 0.6 and Cn = 0.8 compared to results of 
Sec. 4 to orders A^ and A®. The inset shows an blow up of the 
results for small velocities. 



tem of smooth granular particles as long as it remains 
homogeneous. Starting from a pure Gaussian state the 
system develops on a fast time scale to a state where the 
deviations from the Gaussian (described by cumulants) 
are stationary in time and the dynamics is entirely de- 
scribed by a decreasing kinetic energy. 

More technically, we determined formally the full dy- 
namics of the homogeneous cooling state in terms of the 
dynamics of the temperature T(t) and the time depen- 
dent coefScients ai (t) of an expansion of the velocity dis- 
tribution function in Generalized or Associated Laguerre 
polynomials around the Gaussian state. We obtained an 
infinitely large system of non linear ordinary differential 
equations, which can be solved numerically under the as- 
sumption that higher coefficients do not contribute. 

Analytically, we found two main results, i) As far as 
dynamics is concerned, the HCS is characterized by the 
fact that only a few collisions per particle are necessary 
to reach a state where the coefficients are stationary in 
time. Then the entire time dependence is given by a slow 
algebraically decay of the temperature obeying = 
—2ujQ"fT, with 7 depending on all the asymptotic values 
of the coefficients a/, ii) As far as the asymptotic values 
of the coefficients are concerned, the expansion seem to 
converge in the sense that |a;+i| < \ai\ for e„ > 0.6 and 
for this range of e„ we do not find significant differences 
between orders. There exist further stationary but dynam- 
ically instable solutions of the considered differential equa- 
tion, which are far away from the assumption of absolutely 
decreasing and small coefficients, so we cannot make any 
prediction for that cases. For e„ < 0.6 the perturbation 
procedure seems to fail, the assumption we made to trun- 
cate the system of differential equations are severely vi- 
olated and we find a strong dependency on the order of 
approximation. We have no answer if going to higher or- 
der or choosing a more suitable truncation scheme would 
show that the perturbation procedure nevertheless works, 
or, on the other hand, if the expansion presented here is 
only of an asymptotic type. A reasonable conclusion is 
that the system develops to a state which is very far from 
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a Gaussian and might be better described by an expansion 
around an exponential as discussed in and . 

Although much numerical work has been done on the 
HCS and clustering regimes (see, e.g. [Q, j^, 0, |23 and 
p6|), for the first time event-driven simulations are used in 
the present work to investigate deviations from the Gaus- 
sian distribution in the HCS. Mainly, three aspects were 
considered: 1) The dynamics and asymptotic values of the 
coefficients, 2) the influence on the decay rate of the tem- 
perature, 3) the shape of the velocity distribution func- 
tion. 

1. As long as the system remains in the homogeneous 
state the dynamical behavior as well as the static value 
of a2 could be confirmed. Nevertheless, the statistics 
was too poor to distinguish if higher order analytical 
results give better values. 

2. Similarly, the decay of the temperature was measured 
showing deviations with respect to the analytical re- 
sults and opening the possibility to study other effects 
as correlations. Here we had to assure that the sys- 
tem remains in the homogeneous state, restricting the 
number of particles and therefore the quality of the 
statistics. 

3. Measuring the full velocity distribution function we 
found that the distribution function can be described 
very well by the expansion around the Gaussian as long 
as e„ > 0.6. For smaller e„ the high-energy tails show 
an exponential shape confirming previous results found 
by analytical theory and DSMC simulations . 

A possible extension of our work are systems of rough 
spheres with constant coefficient of restitution or Coulomb 
friction. Strong deviations from the Gaussian are observed 
in the angular velocity distribution function, surprisin gly 
for the cases where the particles are almost smooth |p7| , 

. It would be interesting to perform a similar dynamical 
analysis along these lines. 
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